Modeling Elasticity in Crystal Growth 
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A new model of crystal growth is presented that describes the phenomena on atomic length and 
diffusive time scales. The former incorporates elastic and plastic deformation in a natural manner, 
and the latter enables access to times scales much larger than conventional atomic methods. The 
model is shown to be consistent with the predictions of Read and Shockley for grain boundary 
energy, and Matthews and Blakeslee for misfit dislocations in epitaxial growth. 
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The appearance and growth of crystal phases occurs in 
many technologically important processes including epi- 
taxial growth and zone refinement. While a plethora of 
models have been constructed to examine various aspects 
of these phenomena it has proved difficult to develop a 
computationally efficient model that can be used for a 
wide range of applications. For example, standard molec- 
ular dynamics simulations include the necessary physics 
but are limited by atomic sizes (A) and phonon time 
scales (ps). Conversely continuum field theories can ac- 
cess longer length (i.e., correlation length) and time (i.e., 
diffusive) scales, but are difficult to incorporate with the 
appropriate physics. In this paper a new model is pre- 
sented that includes the essential physics and is not lim- 
ited by atomic time scales. 

To illustrate the features that must be incorporated it 
is useful to consider two examples. First consider the nu- 
cleation and growth of crystals from a pure supercooled 
liquid or vapor phase. In such a process small crystallites 
nucleate (heterogeneously or homogeneously) and grow 
in arbitrary locations and orientations. Eventually the 
crystallites impinge on one another and grain boundaries 
are formed. Further growth is then determined by the 
evolution of grain boundaries. Now consider the growth 
of a thin crystal film on a substrate of a different crystal 
structure. Typically the substrate stresses the over-lying 
film which can destabilize the growing film and cause an 
elastic defect-free morphological deformation [0j|], plas- 
tic deformation involving misfit dislocations 0] , or a com- 
bination of both. Thus the model must be able to nu- 
cleate crystallites at arbitrary locations and orientations 
and contain elastic and plastic deformations. While all 
these features are naturally incorporated in atomistic de- 
scriptions, they are much more difficult to include in con- 
tinuum or phase field models. 

Historically, many continuum models have been devel- 
oped to describe certain aspects of crystal growth and 
liquid/solid transitions in general. At the simplest level, 
'model A' in the Halperin and Hohenberg B| classifica- 
tion scheme has been used to describe liquid/solid transi- 
tions. This model treats all solids equivalently and does 
not introduce any crystal anisotropy. Extensions to this 
basic model have been developed to incorporate a solid 



phase that has multiple states that represent multiple 
orientations Kpl or recently M an infinite number of ori- 
entations. Unfortunately these models do not properly 
include plastic and elastic deformations. Other models 
p-[l2| have sought to include elastic and plastic defor- 
mations with reference to a specific reference lattice, but 
cannot account for multiple orientations. 

In this work these limitation are overcome by con- 
sidering a free energy that is minimized by a periodic 
hexagonal (i.e., solid) state. Such free energies have 
arisen in many other physical systems [^,0 (such as 
water/surfactant systems, copolymers, Rayleigh-Benard 
convection and ferromagntic films) and in some instances 
have even been described in terms of crystalline termi- 
nology |14|. The model used in this work describes the 
statics and dynamics of a conserved field, f/;, by the fol- 
lowing free energy and equation of motion. 
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dip/dt = V^ (ST/Sil;) + T], 
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where 77 is a stochastic noise, with zero mean and corre- 
lations < T]{r,t)r]{r',t') >= -GV^5{f- r')S{t - t'), and 
G — hereafter, Qo and e are constants. The field -0 rep- 
resents the average atomic positions. The focus of this 
paper is on two dimensions (2D); it is straightforward to 
extend these calculations to 3D. In two dimensions, this 
free energy is minimized by striped (i/'s), hexagonal (iph) 
and constant (ipc) states depending on the average value, 
Tp, oi if). To estimate the phase diagram these states can 
be approximated f/'s = A^ su\{qax) + ^ 
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and ipc — Tp- Substituting these expressions into Eq. 
(|l|) and minimizing subject to the constraint J dfijj — ip 
gives the values for the characteristic amplitudes A and 



wavenumber q |]i 



and the phase diagram show in Fig. 
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FIG. 1. a) Mean field phase diagram. In this figure the 
'hatched' regions are coexistence regions, b) Grain boundary 
energy. The points are from numerical simulations and the 
line a fit to the Read-Shockley equation, c) Grain growth. In 
this figure the number defects is plotted and the solid line is 
a guide to the eye. d) Epitaxial growth. The points are from 
numerical simulations and the lines are best fits. 

The linear elastic properties of the hexagonal phase 
can be determined by calculating the change of energy 
under shear, bulk or dilational deformations within the 
two mode approximation given in Eq. (0) . It is straight- 
forward to find the elastic moduli for the isotropic solid: 
Ci2 = C44 = Cii/3, where 



Ci 



3^-+ V15e-36V'2 qt 



/75. 



(4) 



For these coefficients [n8[ the Poisson ratio is i^ = 1/3 
and the shear modulus is fj, = C12. 

The energy per unit surface length, E^, between grains 
that differ in orientation by an angle 9 was determined 
numerically and, in Fig. (|l|b), compared with the pre- 
diction of Read and Shockley (l^, i.e., El = EmO[1 — 
lii{9/9M)] where Em and 9m are constants. The parame- 
ters of the simulations were [t.'ip.qo) = (4/15, 1/4, 1). In 
all the simulations conducted in the paper the time and 
space size were At = 0.01 and Aa; = 7r/4 respectively. 
The Read-Shockley form closely fits the data for 9m — 
27.85° and Em = 0.064. The maximum angle, 9m is sim- 
ilar to those observed in experiment [|l7| . The maximum 
energy /length, Em, can be estimated using the Read- 
Shockley equation |1^ in 2D [Q using the elastic con- 
stants estimated above, Em = /^(l + i')b/A-K = Ci2&/37r, 
where h = 2tt /q^ is the magnitude of the Burger's vec- 



tor. For the parameters used in the simulations this gives 
Em — 0.044, as compared to our measured value of 0.064. 
The main advantage of the current approach over 
molecular dynamics simulations is the time scales that 
are accessible. To illustrate this point it is useful to 
calculate the diffusion time using a standard Block- 
Floquet linear analysis. In such an approach the dy- 
namics of a perturbation (^"0) around an equilibrium 
crystal state {ipeq) is obtained by first substituting ^ = 
Stp + tpeq into Eq. (^) and linearizing in Sip. Substi- 
tuting appropriate forms for the equilibrium state (i.e., 
ipeq — ip + J2 ^n,me'xv[iinqxX + mqyy)]) and perturba- 
tion along one of the three principle axes (i.e., Sip = 
Yl ^o,n,mit)exp[i{n{qx+Q)x+mqyy]) and integrating over 
exp[i{kqxX + Iqyv)] gives an equation of motion for the 
modes San^m- Using the approximation for tpeq given 
in Eq. (1^) and four modes (i.e., {Sai^i, Sa^i^i, 6ai^-i, 
(5a_i__i) to describe the perturbation leads to a set of 
four linear ordinary differential equations that can be 
solved analytically, assuming San^m ^ exp(ciji). The 
largest eigenvalue is equal to u — Sq^Q'^ implying a 
diffusion constant of Sq^. In terms of times steps, this 
implies it takes roughly 1000 times steps for a diffusion 
time, Tjj p6| , for the simulation parameters used in this 
paper. 
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FIG. 2. Vacancy diffusion. In this plot the gray scale is 
proportional to the difference between a perfect lattice and 
the one containing the vacancy at times I/td = 0, 5, 10 and 
15 in a), b), c), and d) respectively. For clarity the magnitude 
of the gray scale was maximized in each figure. 

To confirm these calculations, vacancy diffusion was 
studied numerically by removing a single 'particle' from 
a perfect hexagonal lattice, for the parameters used in the 
grain boundary simulations. The dynamics of a vacancy 
is depicted in Fig. (g). It is important to note that this 
model describes long times scales so the vacancy does not 



'jump' from site to site, but the probability of finding the 
vacancy diffuses into the background matrix as shown in 
Fig. ^ . The rate of spread can be used to numericaUy 
determine the diffusion constant and was D — 2.96 which 
is quite close to the approximate calculation (i.e., D — S, 
for Qo = 1). 




FIG. 3. Grain Growth. Gray scale plot of the order param- 
eter ^ at t/m = 700. This plot contains one sixty-fourth of 
the full simulations cell. 

The ability of the continuum model to describe multi- 
ple crystals in arbitrary orientations and locations with 
the appropriate grain boundaries energies on diffusive 
time scales makes it ideal for the study of grain bound- 
ary growth. To study this phenomena a system of 
size 4096Aa; x 4096Aa; was simulated with parameters 
(e, "0, go) = (0.1,0.05,1). To begin the simulations, 256 
crystals were nucleated at arbitrary locations by placing 
large fluctuations in ip at each nucleation site. As time 
evolves the small crystallites grow from the initial seeds 
until impingement. Eventually the entire systems is filled 
with the hexagonal phase and further evolution contin- 
ues by motion at grain boundaries. At this stage in the 
simulations there were approximately 236,000 particles. 
To study the subsequent dynamics the number defects, 
where a defect is defined as a particle that does not have 
six nearest neighbors, was monitored. The number of de- 
fects is shown in Fig. (np) as a function of the logarithm 
of time. Initially (log(t) < 1.7) the number of defect in- 
creases as the total droplet surface area increases. When 
the droplets impinge, the number of defects begins to 
decreases as local rearrangements take place. At later 
times (log(i/T£)) > 2.0) the grain boundaries evolve at a 
very slow rate. In these simulations it is found that the 
number of defects decreases logarithmically at late times. 

The inclusion of elastic and plastic deformations al- 
low the study of morphological instabilities in epitaxial 



growth. When a film grows on a bulk material that has a 
different lattice constant the film can become corrugated 
("buckle") in an attempt to relieve the strain [nl. The 
buckling of the surface relieves the strain in some regions 
but increases the strain in others. At these locations dis- 
locations eventually nucleate. The critical height. He, at 
which these nucleate is well described by the Matthews 
and Blakeslee equation ||3|| which has the functional form. 
He — Ho[l + ln{Hc/b)]/£, where Ho is a constant and £ 
is the mismatch between the film and bulk lattice param- 
eters, i.e., £ = {afiim - abuik)/abuik, where a is a lattice 
constant. 
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FIG. 4. Epitaxial Growth. In this figure the film and bulk 
particles (defined as local maxima of tp) are plotted as open 
and closed circles respectively. Defects are plotted as solid 
squares. 

To study this phenomena a thin film with qo = qj 
was grown on a bulk sample {qo — 1) with the param- 
eters (£,-0) = (1/4,1/4) for various values of qj. The 
system had a width of 8192A2; and height 1024Aa;. A 
small portion of a sample simulation is shown in Fig. 
(^ for a lattice mismatch of 6.4% , (i.e., £ — {^n/qf — 
2tt / qo) / (2Tr / qo) — 0.064). The buckling phenomena is 
shown in Fig. (Qa) and the nucleation of dislocations in 
Fig. (Hd) . The critical height was defined as the height at 
which the interface velocity changes (i.e., since the num- 
ber of particles arriving at the surface is conserved the 
velocity of the film interface changes when dislocations 
appear and change the total number of interface parti- 
cles). The numerical results for He shown in Fig. (|l]d) 
are consistent with the functional relationship proposed 
by Matthews and Blakeslee ||^. 

The one component model described by Eqs. (M and 
(g) does not support other metastable crystal phases and 
thus cannot be used to, for example, study structural 
phase transitions. However it is straightforward to ex- 



tend the model to include more than one kind of particle 
which can give rise to other metastable crystal phases. 
For example a binary system can be easily modeled by a 
free energy of the form; 



T 



dr{yJi M + Vf - ei] Vi/2 + Vi/4 



+ ^2 [{ql + V')2 - €2] ^-2/2 + V4/4 + ai/^iV'2) (5) 



and the equations of motion; dipi/dt = FiV^iJJ-'/J^/'i+r/i 
and dip2/dt — T2^'^5T /5^2 + V2, where a is a couphng 
constant. The properties (i.e., lattice constant, bulk com- 
pressiblity, etc.) of the individual atoms are controlled 
by the parameters with subscripts 1 or 2 and by the aver- 
age value of ipi and '02 ■ When an individual binary alloy 
droplet is grown a hexagonal pattern typically emerges. 
However when a random initial condition is used the pat- 
terns usually contain more than one crystal phase since 
the energy of the hexagonal state is very close to a face 
centered cubic. One such a configuration is shown in Fig. 
(H). Eventually the system will evolve to a hexagonal 
state. Thus the model can be used to study structural 
phase transitions. 
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FIG. 5. Binary Alloy. In this figure the open circles and 
solid squares are the maxima of tpi and -02 respectively. Two 
different crystal structures have been highlighted by joining 
nearest neighbors of the same species. 

In conclusion, the preceding simulations and calcula- 
tions have provided ample evidence that the model de- 
scribed by Eqs. (|^) and (|^) provide an adequate descrip- 
tion of crystal behavior on long time and atomic length 
scales. Thus this model should provide a useful method 
of studying the crystallization of pure and multicompo- 
nent materials. 
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